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Abstract 

In lattice QCD it is possible, in principle, to determine the parameters in the effective chiral 
lagrangian (including weak interaction couplings) by performing numerical simulations in 
the e-regime, i.e. at quark masses where the physical extent of the lattice is much smaller 
than the Compton wave length of the pion. The use of a formulation of the lattice theory 
that preserves chiral symmetry is attractive in this context, but the numerical implemen- 
tation of any such approach requires special care in this kinematical situation due to the 
presence of some very low eigenvalues of the Dirac operator. We discuss a set of techniques 
(low-mode preconditioning and adapted-precision algorithms in particular) that make such 
computations numerically safe and more efficient by a large factor. 



1. Introduction 

At low energies the physics of the light pseudo-scalar mesons can be described by an 
effective chiral cr-model with SU(3) left x SU(3) right symmetry group. Purely strong 
interaction phenomena and electroweak transitions (such as the non-leptonic kaon 
decays) are both covered by the effective theory if the appropriate interaction terms 
are included in the lagrangian. The challenge is then to compute the associated 
coupling constants from the underlying field theory of the strong and the electroweak 
interactions. 

In the effective chiral theory the physical amplitudes are obtained in the form of 
an asymptotic expansion in powers of the quark masses and the meson momenta. 
The coupling constants appear in the coefficients of this expansion, and to determine 
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their values it is therefore essential to choose a computational strategy where the 
limit of small masses and momenta can be safely reached. 

The so-called e-regime of QCD [1-6], combined with a formulation of lattice QCD 
in which chiral symmetry is exactly preserved [7-15], provides a framework that may 
prove to be particularly suitable in this context. Some encouraging results have in 
fact already been obtained along these lines in the case of the quark condensate [16- 
19] which is one of the free parameters in the chiral lagrangian. In general numerical 
simulations in the e-regime are technically demanding, however, because the lattice 
Dirac operator tends to be ill-conditioned. The problem is linked to the spontaneous 
breaking of chiral symmetry [6,20-27] and is hence present independently of how 
precisely the lattice theory is defined. 

In this paper we discuss a set of fast numerical methods for the computation of the 
quark propagator and the zero-modes of the Dirac operator. Being efficient is clearly 
very important here, but we wish to emphasize that an ill-conditioned system also 
gives rise to questions of numerical stability and accuracy [28] that must be properly 
dealt with. 

Although some of our techniques are more widely applicable, we focus on the case 
of Neuberger's lattice Dirac operator [13] in this paper. After a brief review of the 
e-regime in sect. 2, we then first discuss a particular numerical implementation of 
this operator, following previous work on the subject [15,29-32], which is uniformly 
accurate to a specified level of precision (sects. 3-5). The ability to rigorously con- 
trol the approximation error is seen to be very useful in sect. 6, where we describe 
an efficient way to calculate the index of the Dirac operator, and it also provides the 
basis for the adapted-precision inversion algorithm discussed in sect. 9. Before this, 
in sects. 7 and 8, we show how the zero- mode contribution to the quark propagator 
can be safely separated and introduce a method, referred to as low-mode precon- 
ditioning, that takes care (to some extent) of the potentially very large condition 
number of the Dirac operator in the subspace orthogonal to the zero-modes. 



2. The e regime of QCD 

In numerical simulations of lattice QCD the adjustable parameters are the gauge 
coupling, the quark masses and the lattice size. Usually a euclidean T x L 3 lattice 
is considered, with periodic boundary conditions, where T and L refer to the time 
and the space directions respectively. We shall also assume in this paper that the 
lattice Dirac operator satisfies the Ginsparg- Wilson relation so that chiral symmetry 
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is exactly preserved on the lattice [14]. The e-regime is then reached at fixed gauge 
coupling by scaling T and L to values much larger than the confinement radius and 
simultaneously the light quark masses m = m u ,md,m s to zero in such a way that 
the combination 



remains fixed. In this equation the parameter 

£ = lim lim \(uu}\ (2.2) 

m— >0 V^oo 

is the (unrenormalized) u-quark condensate in infinite volume.f 

2.1 Effective chiral theory 

To leading order the purely strong interaction part of the euclidean action density 
of the effective theory reads 



where U denotes the SU(3)-valued chiral field and M = diag (m u , m^, m s ) the quark 
mass matrix. The normalizations chosen here are such that the parameters F and 
£ coincide with the pion decay constant and the chiral condensate at tree-level of 
the chiral perturbation expansion. 

In the chiral limit the effective theory is expected to provide a correct description 
of the low-energy properties of lattice QCD, up to lattice spacing effects, if F and 
£ are set to the proper values in units of the lattice spacing. These parameters can 
thus be determined by calculating any suitable correlation functions on the lattice 
at small masses and momenta and by comparing the results with the predictions 
of the effective theory. In particular, such computations may be performed in the 
e-regime, where the quark mass term breaks chiral symmetry only very weakly and 
where the zero-momentum modes U{x) = V dominate the partition function of the 
effective theory [1]. There are many different ways to study the matching between 
the effective theory and lattice QCD in this regime, most of which still need to be 
explored. 

f For simplicity we do not consider the case of quenched or partially quenched QCD in this section. 
Moreover we stick to a continuum notation and drop all factors that formally converge to 1 in the 
continuum limit. 



X = m£V, 



V = TL 3 



(2.1) 



C = \F 2 tr {d^d^U} - \Y.tT{UM + MW} , 



(2.3) 
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An example that we wish to quote for illustration is the two-point function of the 
left-handed current 

31 = W^P-^P, P± = W± 75), (2.4) 

which is represented by 

j« = \F 2 tr{\ a Ud^} (2.5) 

in the effective theory. Standard notations are being used here, where ip denotes 
the quark field with flavour components (u,d, s, . . .) and A a an SU(3) generator that 
acts on the flavour indices of the fields. At non-zero time separations and vanishing 
spatial momentum, this correlation function turns out to be proportional to F 2 /T 
[4,5,33] and it is thus expected to be particularly suitable for the calculation of the 
parameter F. 

2.2 Weak-interaction couplings 

The effective couplings associated with hadronic weak transitions can perhaps be 
computed in a similar way by matching correlation functions in the e-regime. In the 
case of the CP conserving K — > tttt decays, for example, the weak interaction density 
that is obtained after integrating out the VF-bosons is, to a good approximation, 
equal to a linear combination of [34-38] 

± = {(s ltl P_u)(u^P.d) ± (s^P.d)(u lfl P.u)} - (u - c) , (2.6) 

O m = {m 2 c - m 2 u ) {m d {sP + d) + m s {sP_d)} . (2.7) 

The operator O m is usually omitted from this list, because it does not contribute in 
transition matrix elements where the total energy-momentum is conserved. Mixings 
of and O m with other operators are, incidentally, excluded by the exact sym- 
metries of the lattice theory. 

In the effective chiral theory, the operators and O m are expected to be repre- 
sented by certain linear combinations of the fields [39] 

01 = lF 4 {(Ud^) 13 (Ud^) 21 + 

(Ud^hsiud^u + l^Ud^hs}, (2.8) 

02 = lF 4 (d^Ud^) 23 , (2.9) 



4 



3 = iF 2 £([/M + Mt[/t) 23 , 



(2.10) 



and our task is then to determine the dimensionless coefficients in these linear com- 
binations. An obvious possibility in the e-regime is to match three-point correlation 
functions such as 

C ab (x ,yo) = / L d 3 xd 3 y<j;(x)O ± (0)^(y)> (2.11) 
Jo 

plus a similar correlation function with only one left-handed current (so as to be able 
to disentangle the contributions of the octet operators O2 and O3). Explicit com- 
putations in chiral perturbation theory suggest that these correlation functions are 
indeed well-behaved and that, in principle, the desired coefficients can be calculated 
along these lines [40]. 

2. 3 Spectrum of the Dime operator 

The last ten years or so have seen a remarkable development, involving both random 
matrix and chiral perturbation theory, that culminated in the analytical determi- 
nation of the distributions of the low-lying eigenvalues of the Dirac operator in the 
e-regime [6,20-27] (for a review see ref. [41]). An important insight gained in the 
course of this work is that the eigenvalues scale like (SF) _1 at large volumes. In 
particular, the spectral gap in the subspace orthogonal to the chiral zero-modes is, 
on average, roughly of this size. 

For the calculation of the quark propagator in the e-regime, the presence of some 
very low eigenvalues (an effect that has been confirmed numerically on small lattices 
[42,43,19]) is a source of difficulty. To make this a bit more concrete, let us consider 
a 32 4 lattice with spacing a = 0.1 fm and let us assume that S is around (250 MeV) 3 
on this lattice. From the discussion above we then deduce that the expectation value 
of the spectral gap is about 5 x 10~ 4 in lattice units. Note that the quark masses are 
scaled proportionally to (EV) _1 in the e-regime and thus do not provide a strong 
infrared cutoff on the spectrum. 

The computation of the quark propagator, using the conjugate gradient algorithm 
for example, is perhaps still possible under these conditions, but may require a very 
large number of iterations. Moreover, the accuracy of the solution vector that can be 
reached on any given machine may be unsatisfactory, especially if the propagator is 
calculated via the so-called normal equations (as is the case if the conjugate gradient 
algorithm is used) [28]. Effectively these technical problems set a limit on the space- 
time volumes where numerical simulations can be performed and thus on how close 
to the chiral point one can get in practice. 
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3. Lattice Dirac operator 



We set up the lattice theory as usual on a finite four-dimensional lattice with spacing 
a and periodic boundary conditions. Our notational conventions are summarized in 
appendix A. Following refs. [7-14] we consider a formulation of lattice QCD where 
the lattice Dirac operator D satisfies the Ginsparg- Wilson relation 

75D + ,075 = aD~f 5 D, = 75D75, (3.1) 

for some positive constant a proportional to a. In the presence of a given background 
gauge field, the propagator of a quark with mass m is then given by 

<^(x)?(y)> = {Z)- 1 }(x,y), (3.2) 

where {. . -}(x, y) stands for the kernel in position space of the operator in the curly 
bracket and 

D m = (1 - \am)D + m (3.3) 

is the massive lattice Dirac operator. We shall always assume that < am < 2. 

As already mentioned in sect. 1, we concentrate on the case of the Neuberger-Dirac 
operator in this paper, even though many of the techniques that will be described 
are more generally applicable. Explicitly this operator is given by [13] 

^ = ^{l + 7ssign(Q)}, (3.4) 

Q = 75 (aD w - 1 - s) , a = — , (3.5) 

1 ~r S 

where s is an adjustable parameter in the range |s| < 1 and D w denotes the standard 
Wilson-Dirac operator (appendix A). 



4. Minmax polynomial approximation 

For the numerical implementation of the Neuberger-Dirac operator we need to find 
an algorithm that evaluates sign(Q)?7 on any given fermion field r] to a specified 
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precision. In the present section we take a first step in this direction by constructing 
an optimal polynomial approximation to sign(x) in a range that excludes a small 
interval around the origin. More precisely we are looking for a polynomial P(y) of 
degree n that minimizes the error 

5 = max |%)| , h(y) = 1 - y/yP(y), (4.1) 

e<y<l 

for specified e > 0. In the range \fe < \x\ < 1 the function xP{x 2 ) then approximates 
sign(x) uniformly with a maximal deviation 5. 

4-1 Existence and uniqueness 

Polynomial approximations that minimize the maximal relative error are referred 
to as minmax polynomials. They can be shown to exist and to be uniquely deter- 
mined provided the function that is to be approximated does not vanish. Moreover, 
practical methods have been devised, based on the characteristic properties of the 
associated error function, that allow one to compute them reliably (see ref. [44], for 
example, or any other good book on approximation theory). 

For reasons of numerical stability, it is advantageous to represent the polynomial 
in the form of a series 

n 

P(y) = c kT k (z), z = (2y-l- e)/(l - e), (4.2) 

k=0 

of Chebyshev polynomials (appendix A). Our task is then to adjust the coefficients 
Cfc so that the deviation (4.1) is minimized. 

4-2 Error bounds 

The characteristic property alluded to above of the minmax polynomial is that the 
error function h(y) has (at least) n + 2 extrema of equal height and alternating sign, 
i.e. h(y) is oscillating around zero with constant amplitude (fig. 1). Evidently it is 
not possible to know in advance where these extrema are located. 
Let us now choose an arbitrary sequence 

e < yo < yi < ■ ■ ■ < y n +i < i (4.3) 

of n + 2 values. We can then find a polynomial such that 

h(yi) = (-l) l u for all Z = 0,...,n + 1. (4.4) 



7 



_2 I I I I I I I I I I I 

-1 -0.5 0.5 1 

Fig. 1. Minmax polynomial approximation of sign(x) for e = 0.0025 and n — 22 
(full line). In the range y/e < \x\ < 1 the relative error of this approximation is 5%, 
while for an error of 0.1% a polynomial of degree n — 57 is required (dashed line). 

In fact these are just n + 2 linear equations in the unknowns Co, . . . , c n , u. This poly- 
nomial has the alternating property on the chosen set of points, but not necessarily 
so on the approximation interval. It can be shown, however, that 

\u\ < max |/i*(w)| < max 
e<2/<! «<y<i 

where h*(y) denotes the error function of the minmax polynomial. The solution of 
the linear equations (4.4) thus yields a nearly optimal polynomial if these bounds 
are tight. 

4-3 Maehly's exchange algorithm 

The basic idea of this algorithm is to start from some well-chosen initial set of points 
(4.3) and to adjust them subsequently so as to tighten the bounds (4.5). This can be 
achieved, for example, by choosing the new points to be the nearby extremal points 
of the current error function. 

A good initial set of points for this recursion are the extremal points of the Cheby- 
shev polynomial T n+1 (z). Only few iterations are then required until a polynomial 
is obtained that is very nearly optimal. There is in fact little to be gained by per- 
forming further iterations once the lower and the upper bound in (4.5) differ by no 
more than (say) 0.01 |u|. Evidently this procedure is extremely robust, because the 



(4.5) 



quality of the current polynomial (with respect to the best possible approximation) 
can be rigorously controlled. 



4-4 Implementation details 

Inspection shows that the linear system (4.4) is well-conditioned and it can thus be 
solved straightforwardly by QR factorization, using Householder transformations, 
followed by backward substitution (see ref. [28] for example). 

To localize the extrema of h(y), a few bisection steps are sufficient in general, 
since there is no need to be very precise at this point. An important detail here is 
that many digits are lost when calculating the error function, even if the Clenshaw 
recursion is used to evaluate the Chebyshev series (4.2) [45]. This becomes a problem 
when bjYln approaches the machine precision, and with standard 64 bit arithmetic 
the approximation accuracy 5 that can be reached is thus limited to values greater 
than n x 10~ 15 or so. 

4-5 Efficiency of the approximation 

The precision of the approximation to the sign function provided by the minmax 
polynomial scales approximately like 



where A = 0.41 and b = 2.1. At large values of roy^, this empirical formula is fairly 
accurate, but it slightly over-estimates the actual value of 6 when the degree of the 
polynomial is less than 30 or so. With respect to the Chebyshev approximation that 
was used in ref. [15], the precision is improved by about a factor of 2. 

We finally mention that, on a current PC processor and for degrees n < 400, only 
a few seconds are required to obtain the coefficients Ck of the minmax polynomial. 
It is hence practical to recompute these polynomials whenever they are needed. 



5. Uniform approximation of the Neuberger Dirac operator 

The minmax polynomial P n ,e(y) of degree n on the interval [e, 1] yields an approxi- 
mation of the operator sign(Q) through 



5 = Ae~ bn ^ 



(4.6) 



sign(Q) ~XP„, e (A 2 ), X = Q/\\Q\\. 



(5.1) 
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If e is chosen so that Q 2 > e||Q|| 2 , the error in this formula is an operator with norm 
less than or equal to S. In other words, the approximation error is always bounded 
by 5\\rj\\, uniformly in the field r\ to which the operator is applied. 

This method is straightforward but actually not recommendable if Q 2 has some 
exceptionally low eigenvalues (as is often the case in practice [15]). It is far more effi- 
cient under these conditions to first separate the few lowest modes and to treat them 
exactly. Evidently this should be done in such a way that the total approximation 
error remains under control. 

5.1 Low-mode projectors 

The spectrum of Q in the vicinity of the origin can be reliably determined by mini- 
mizing the Ritz functional of Q 2 [46,47]. This technique also yields an approximation 
to the associated eigenvectors, and we now need to discuss by how much these vectors 
deviate from the true eigenvectors. 

So let us assume that a specified number / of approximate eigenvectors has been 
computed. We denote by V the linear space spanned by these vectors and by P the 
corresponding orthonormal projector. It is then trivial to calculate the parameter 

q= max ||(l-P)Qv|| (5.2) 

v€V,\\v\\ = l 

which measures the deviation of V from being an exact eigenspace of Q. The im- 
portance of this parameter is clarified by 

Lemma 5.1. Let v\, . . . , v\ be the eigenvalues of FQF in the subspace V . Then 
there are I linearly independent eigenvectors of Q with eigenvalues Ai, . . . , A; such 
that \uk — Afc| < Q for all k = 1, . . . , I. 

The proof of the lemma (which holds independently of how V was obtained) is given 
in appendix B. 

In the following we shall take it for granted that the eigenvalues u k are separated 
from zero and from the rest of the spectrum of Q (all eigenvalues A except Ai, . . . , A;) 
by a distance greater than q. So far we have in fact never encountered a situation 
where this condition could not be satisfied, and we thus omit any further discussion 
of this point. The presence of a spectral gap around zero implies that the subsets of 
positive and negative eigenvalues can be identified without any numerical ambiguity. 
If we introduce the eigenvectors 

u k £ V, FQu k = v k Uk, (u k ,Uj) = 5 k j, (5.3) 
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the associated orthonormal projectors are given by 

p+ = Yl Uk ® K) f and p - = E Uk ® ( Ufc ) f ( 5 - 4 ) 

respectively. In the same way we may define the projectors (P±) eX act to the subspaces 
spanned by the corresponding exact eigenvectors of Q (whose existence is guaranteed 
by lemma 5.1) with positive and negative eigenvalues Afc. 

We may now ask how accurately the computed projectors P± approximate the 
exact projectors (P±) exact- The answer to this question depends on the size of the 
residues 

Qk = \\{Q ~ v k )u k \\ (5.5) 

and also on the distances between the eigenvalues of Q. A small value of g k alone 
does, in fact, not exclude sizeable mixings of u k with several eigenvectors of Q if 
these have eigenvalues that are within a distance g k of v k . 

Eventually we are interested in estimating the deviation of the projectors rather 
than that of the individual eigenvectors, and what counts in this case is the distance 
dk of Vk from the exact spectrum of Q in the subspace that is orthogonal to the 
range of (P + ) exact if > or (P_) exact if v k < 0. The quality of the approximation 
is then controlled by the parameters 

4= E Sl/d 2 kl k±>0, (5.6) 
as the following lemma shows. 

Lemma 5.2. Provided the spectral conditions listed above are satisfied, the bound 
holds where it is assumed that 2(1 + 1)k±(1 + 2k±) < 1. 

The proof of the lemma is given in appendix C. Here we only note that the param- 
eters k± are very small in the cases of interest. The right-hand side of eq. (5.7) is 
then practically equal to k± and the condition on the last line is trivially fulfilled. 
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5. 2 Approximation formula and error bound 

It is now straightforward to write a program that computes the projectors to a 
specified accuracy. In this program the number of low modes to be included in 
the projectors should be determined dynamically in such a way that the spectral 
distance from the other modes is not accidentally very small. The parameters n± 
can then be estimated without difficulty and the minimization of the Ritz functional 
is stopped when the desired level of precision is reached. 

Once the projectors are computed, eq. (5.1) may be replaced by 

sign(Q)~P + -P_ + (l-P+-P_)XP n>e (X 2 ), X = Q/\\Q\\, (5.8) 

where e is set to a value equal to (or perhaps slightly less than) the smallest eigen- 
value of Q 2 /||(3|| 2 in the sector orthogonal to the low modes. The properties of the 
minmax polynomial, lemma 5.2 and simple triangle inequalities then imply that the 
associated approximation D m to the massive Neuberger-Dirac operator satisfies 

II An - Anil < - (1 + 8 - \am) {2{k + + k_) + 5} (5.9) 
up to terms proportional to k±5 and k± (from now on these will be neglected). 

5.3 Miscellaneous remarks 

In the tests that we have performed we never had any difficulty in reaching accu- 
racy levels k± of 10 -8 or even 10 -10 . Moreover, the actual approximation error is 
significantly smaller than suggested by lemma 5.2, by a factor 10 at least. Since the 
calculation of the projectors consumes a relatively small amount of computer time, 
they will normally be obtained to a precision that is estimated to be sufficient, by 
a wide margin, for the applications that one has in mind. The approximation error 
(5.9) is then essentially determined by the degree n of the minmax polynomial. 

A numerically safe method to evaluate the polynomial in eq. (5.8) is provided by 
the Clenshaw recursion [45]. The accumulated rounding errors are then on the or- 
der of n times the machine precision. In particular, if standard 32 bit floating-point 
arithmetic is used (which can be profitable at intermediate stages of the algorithms 
described later), a rounding error of about 5n x 10~ 7 should be added to the theo- 
retical approximation error 5. 

We finally note that the degree of the minmax polynomial will usually be signif- 
icantly larger than the number of the separated low modes. The application of the 
projectors in eq. (5.8) then requires only a small fraction of the total execution time 
and the associated rounding errors can also be safely neglected. 
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6. Calculation of the index of the Dirac operator 

6.1 Preliminaries 

The operator D commutes with 75 and thus leaves the subspaces of fermion 
fields with definite chirality invariant. Using the Ginsparg-Wilson relation, it can 
be shown that the action of the operator in these subspaces is given by the hermitian 
operators 

D ± = P±DP± (6.1) 

up to a proportionality factor equal to 2/a. Moreover, the spectrum of the low-lying 
eigenvalues of D + and D~ is exactly the same, including degeneracies, apart from 
the fact that there can be a surplus of zero- modes in one of the sectors. 

It happens only accidentally (with probability zero strictly speaking) that there 
are zero- modes with both chiralities, and we shall thus exclude this exceptional case 
in the following. The index v of the Dirac operator is then given by 

v = ariQ (6-2) 

where a denotes the chirality of the sector that contains the zero-modes and no the 
number of these modes. 

6.2 Zero-mode counting 

The low-lying eigenvalues of D + and D~ can, in principle, be found by minimizing 
the associated Ritz functionals [46,47]. In terms of computer time such calculations 
tend to be very expensive, however, since each application of the Neuberger-Dirac 
operator to a given fermion field involves from say 50 to several 100 applications of 
Q 2 . Any numerical method that requires a precise computation of the eigenvalues 
should hence be avoided. 

The strategy that we propose is to run the minimization program in both chiral- 
ity sectors simultaneously in a controlled way, where the precision is improved in 
steps by a specified factor. After every iteration the current estimate of the lowest 
eigenvalue in each sector provides a rigorous upper bound on the gap above zero in 
the spectrum of D + and D~ respectively. The program then proceeds to lower the 
larger bound and continues in this manner until the estimated relative error on any 
one of the eigenvalues in the two sectors is determined to be less than 10%. At this 
point we know that the zero-modes (if any) must be in the other chirality sector. 
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To count the zero-modes the minimization program must now be run once more in 
that sector, improving the precision on the calculated eigenvalue in steps as before. 
As soon as the eigenvalue drops below the spectral gap in the other sector (which 
was determined to an accuracy of 10% in the previous calculation), there has to be 
at least one zero-mode, independently of the current level of precision, because the 
Ritz functional always provides an upper bound on the true eigenvalue. We can then 
restart the minimization program, requiring a second eigenvalue to be computed, 
and so on. 

This procedure terminates when an eigenvalue is found that is definitely above 
zero, which will be the case if the relative error is estimated to be less than 10% for 
example. It should be emphasized that in this way we never require any eigenvalues 
to be calculated to high precision and yet the zero-mode counting is rigorously 
correct. 

6.3 Using reduced-precision operators 

In the discussion above we have implicitly assumed that the approximation error 
(5.9) of the numerical representation of the operators can be neglected. An 
important point to note is that the associated error on the eigenvalues is at most 



and the same is true for the gradients of the Ritz functionals of D . We may, 
therefore, start the minimization of the Ritz functionals using a relatively poor 
approximation of the Dirac operator, and then gradually decrease 5 so that the 
error u is always smaller than the current magnitude of the gradients. 

Apart from being in control of the numerical errors, an obvious advantage of this 
procedure is that the degree of the minmax polynomial is not larger than absolutely 
necessary, at all stages of the calculation. The total computational effort is thus min- 
imized. Moreover, a further acceleration may be achieved by using single-precision 
arithmetic in the programs for the Neuberger-Dirac operator as long as 5 is not too 
small (cf. subsect. 5.3). On current PC processors this saves almost a factor of 2 in 
execution time [48]. 

6.4 Error estimation 

In the algorithm described in subsect. 6.2, the error on the calculated eigenvalues 
must be estimated in a reliable way. For any given eigenvalue, the error is a sum 
of two contributions, one coming from the numerical approximation of the lattice 




(6.3) 



a, 
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Dirac operator and the other from the fact that, in each step of the procedure, the 
minimum of the relevant Ritz functional is only found to some specified accuracy. 

A reliable but often poor estimate of the latter error is provided by the gradient 
of the Ritz functional [46,47]. However, for reasons of numerical stability it is in 
any case advisable to compute an additional eigenvalue in each sector, with possibly 
reduced precision [47] , and the Temple inequality may then be used to obtain a more 
accurate estimate of the error on the other eigenvalues [49]. 



7. The quark propagator: general strategies 

The presence of chiral zero-modes (which is the normal case in large volumes) com- 
plicates the calculation of the quark propagator, and it is mainly this issue that we 
wish to address in this section. We assume in the following that the quark mass m 
is strictly positive and that a gauge field configuration is being considered where the 
zero- modes (if any) have chirality a = +1. 

7.1 Negative chirality components 

The computation of the quark propagator amounts to solving the linear equation 
D m ^j = 7] (7.1) 

for a set of source fields rj. We first write the solution in the form 

^ = (DiD^DU (7.2) 

and note that the operator in brackets commutes with 75 . The components of the 
solution vector with negative chirality are thus given by 

P_V> = {DlD^P.Dln. (7.3) 

In this equation the inversion of DmD m takes place in the chirality sector that does 
not contain the zero-modes, and the only potential difficulty is then that there can 
be some very low non-zero eigenvalues of D^D in this subspace (cf. subsect. 2.4). 

We shall come back to this problem in the next two sections, where we discuss the 
inversion of D\nD m in the negative chirality sector. For the time being, we merely 
assume that an efficient algorithm is available which allows us to invert the operator 
in this sector, to a specified precision and at arbitrarily small quark masses. 
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7.2 Positive chirality components 

Once the negative chirality component P_V> has been computed, it can be inserted 
on the right-hand side of the equation 

D m P + ij = n- D m P_iP (7.4) 

for the positive chirality components P+tp. A possible expression for the solution of 
this equation is 

P + 4> = (P +J D m P + )" 1 {P +V - P + D m P_^} (7.5) 

which involves an inversion of P + D m P + in the sector that contains the zero-modes. 
It is now very important to note that 

P+D m P + = P + { \a (1 - \am) D^D + m} P + (7.6) 

is a positive operator in this sector whose condition number is at most 2 /am and 
thus comparatively small. If the quark mass in lattice units is not too close to zero, 
the inversion in eq. (7.5) is hence well-conditioned and the solution can be computed 
straightforwardly, using the conjugate gradient algorithm for example. 

7.3 Separation of the zero-modes 

Arbitrarily small quark masses can also be handled but require a subtraction of the 
zero-mode component from the source field. Equation (7.5) then becomes 

= ^PoP+V + (P+D^+r 1 {(1 - P )P +V - P + D m P_4>} , (7.7) 

where Po denotes the projector to the subspace spanned by the zero- modes. Note 
that the field in the curly bracket is orthogonal to this subspace, and the operator 
inversion in eq. (7.7) thus remains well-defined even in the limit of a vanishing quark 
mass. The situation is then practically the same as in the negative chirality sector 
discussed in subsect. 7.1. In particular, the acceleration techniques described in the 
next two sections can be applied here too. 

At this point we still need to explain how to compute the zero-modes (and thus 
the projector Pq). First recall that a basis of approximate zero-modes is obtained 
in the course of the calculation of the index of the Dirac operator (cf. subsect. 6.2). 
This basis may not be very accurate, but we can now easily improve on the accuracy 
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by noting that 



PoX = (1 - \aD) P +X - DiptD^P-rfP+x 



(7.8) 



for any field %. The operator inversion in this formula is in the sector with negative 
chirality (where there are no zero-modes) which we have assumed to be feasible. 

It is our experience that this procedure is numerically stable and yields accurate 
results when applied to the approximate zero- modes. More precisely, the error on 
the calculated zero-modes [after the application of eq. (7.8)] is estimated to be at 
most r + 2oj / 1 A 1 1 in this case, where oj is the accuracy of the approximation to the 
Neuberger-Dirac operator used in the numerator, Ai the first non-zero eigenvalue of 
D and r the relative residuum associated with the numerical inversion of D.\ 

7.4 Left-handed propagator 

The correlation functions that we have proposed in sect. 2 to probe the weak-inter- 
action operators involve only left-handed quark and antiquark fields. To evaluate 
these correlation functions, it thus suffices to compute the purely left-handed com- 
ponents of the quark propagator, which are given by 



From this expression it is immediate that the complications arising from the zero- 
modes are absent in this case. In fact we only need to invert DjnD m in the chirality 
sector where there are no zero- modes. Note that if the zero- modes had the opposite 
chirality, we would write the propagator in the form 



so that the inversion of DlnDm is again in the "good" sector. 

Such correlation functions thus appear to be particularly attractive from the nu- 
merical point of view. Moreover, their behaviour in the chiral limit tends to be less 
singular than is generally the case [50,33,40]. 



I An argumentation similar to the one at the beginning of appendix C shows that the deviation 
|(1 — Po)xll °f anv field x from being an exact zero-mode is less than or equal to ||Z?x||/|Ai |. The 
approximation error can thus be rigorously controlled. 



{P.(DlD m )~ 1 P_DlP + } (x,y). 



(7.9) 



{P.DlP + (DlD m )- 1 P + } (x,y) 



(7.10) 
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8. Low-mode preconditioning 

We now proceed to discuss the linear system 

Ai/j = v , A = DlD m , (8.1) 

in the negative chirality sector, assuming as before that the zero-modes (if any) are 
in the other sector. 

As already emphasized in sect. 2, the presence of some extremely low eigenvalues 
of A forbids a straightforward application of the established inversion algorithms 
to solve the equation. To reduce the condition number of the system, an obvious 
possibility is to calculate the few lowest eigenvalues of A and to compute the com- 
ponents of the solution vector along the associated eigenspace and its orthogonal 
complement separately. However, on large lattices this proposition is not practical, 
because the computation of the eigenvectors to the required level of precision is far 
too expensive in terms of computer time. 

In this section we describe a version of low-mode preconditioning that works out 
even if the eigenvectors cannot be calculated very accurately. Effectively the algo- 
rithm compensates for the approximation errors through a simple block diagona- 
lization that can be implemented exactly. 

8.1 Approximate eigenvectors 

So let us assume that ei, . . . , e n is a set of orthonormal Dirac fields with negative 
chirality that satisfy 

Ae k = a k e k + r k , (e h r k ) = 0, \\r k \\ < uj k a k , (8.2) 

for all k,l and some positive numbers a k and u k . Evidently, if the bounds u; k on 
the residues r k are small, these fields can be regarded as approximate eigenvectors 
of A with eigenvalues a k . In the following we shall also take it for granted that 

(v , Av) > 7 \\v || 2 , 7 = maxafc, (8-3) 

k 

for all negative chirality vectors v in the orthogonal complement of the vector space 
E spanned by e\, . . . , e n . This guarantees that ei, . . . , e n approximate the lowest n 
eigenvectors of A rather than an arbitrary set of eigenvectors. 

It should be emphasized that no further assumptions need to be made and that 
the error bounds u> k , in particular, do not need to be extremely small (in most cases 
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values of 0.1 or so will do). In practice the vectors ei,...,e n can be computed by 
minimization of the Ritz functional of A and stopping the program when the esti- 
mated relative errors of the calculated eigenvalues (as determined from the gradient 
of the Ritz functional [46,47]) drop below the specified values of uJk-] 

8.2 Preconditioned system 

We now introduce the projector P to the subspace E and note that the linear system 
(8.1) is equivalent to 

(1 - P) {A - AP{PAP)~ 1 PA} (1 - P)V> 

= (1 - P)t] - (1 - P)AP(PAPy 1 Pr ] , (8.4) 

PAPip = Pn- PA(1 - P)ip. (8.5) 

The first of these equations determines the component (1 — P)ip of the solution in 
the subspace orthogonal to E. Once it has been calculated, the parallel component, 
Pip, is obtained from the second equation. 

To make this a little more explicit, it is helpful to introduce the field 



<p = (i-P)^ = V-J> fc ( efc >^)- ( 8 - 6 ) 
k=i 

Equation (8.4) may then be written as 

n ^ n 1 

(l-P^-Vr t -(r^) = (l-P)v- y>fc— (e k ,v)- (8-7) 



, . "it r-; o>k 

k=l k=l 



Together with the constraint Pip = 0, this is a well-defined system whose solution 
immediately yields the complete field through 

n 1 

*P = <P + Ve fe — {(e k ,v)-{r k ,<t>)}. (8.8) 

The scalar products in these formulae give rise to a computational overhead that 
is completely negligible compared to the effort required for the application of the 



t Since the eigenvectors of A are mass-independent, we may set m = in this calculation and use 
the same vectors ei, . . . ,e n at all values of m. Only the eigenvalues and the residual vectors rj, 
then need to be recalculated if the quark mass changes. 
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operator A. This assumes, of course, that the decomposition (8.2) has been worked 
out before the conjugate gradient program is started and that it is possible to keep 
the vectors e k and r k in memory during the whole calculation. 

8.3 Condition numbers 

The operator on the left-hand side of eq. (8.7) is the difference of two operators, 



that act in the orthogonal complement of E and that are both positive. The lowest 
eigenvalue of the first operator is greater or equal to 7, while the norm of the 
second is bounded by ^fc =1 (^fc) 2 afc- This implies that the condition number of the 
preconditioned system (8.7) is at most 



which is practically equal to 4/a 2 7 if the error bounds ljj, are smaller than 0.2 x n -1 / 2 
for example. Since the condition number of the original system is about 4/a 2 divided 
by the lowest eigenvalue of A, a reduction in the condition number by roughly a 
factor of maxfc afc/min^ a k is thus achieved. 

Evidently this ratio depends on the gauge field configuration and on the number 
of modes that are being subtracted. With n as small as 4 we have found factors of 
30 and more. An interesting remark in this context is that the improvement factor 
is roughly independent of the lattice parameters, because the low-lying eigenvalues 
of A scale proportionally to (SF) -2 (cf. subsect. 2.3). The probability distribution 
of the ratio of the nth to the first eigenvalue is in fact universal in random matrix 
theory and analytically computable. 

8.4 The low-mode contribution — a precision issue 

The program that minimizes the Ritz functional of A initially yields a set of approx- 
imate eigenvectors e~\ , . . . , e n that are orthonormal to machine precision but that are 
not guaranteed to satisfy eq. (8.2) very accurately. We can, however, easily correct 
for this deficit by defining the vectors e k to be the exact eigenvectors of A in the 
subspace spanned by the given vectors. The numbers a k are then the eigenvalues of 
the matrix 




(8.9) 



(4/a 2 )x{ 7 -ELiK) 2 «4 



-1 



(8.10) 



M M = (e k ,Ae~i) 



(8.11) 
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and the residues are obtained through = (1 — P)Aek- 

While the subspace spanned by the eigenvectors is by definition given to machine 
precision, the calculation of the eigenvalues is a possible source of numerical 
inaccuracy. Note that the error in these numbers is amplified in eq. (8.8) and may 
thus easily affect the precision of the solution vector. 

The matrix Mki should therefore be calculated using a high-precision approxima- 
tion to the Neuberger-Dirac operator. It is also important to apply the operator A 
in its original form, eq. (8.1), so that the eigenvalues a?c of the matrix are obtained 
with an error equal to the approximation error on D m times (o^) 1 / 2 (rather than 
just the approximation error). 



9. Adapted-precision inversion algorithm 



Iterative improvement combined with the use of low-precision arithmetic is men- 
tioned in many books on numerical mathematics as an effective method to speed 
up the computation of the solution of linear systems (see refs. [28,45] for example). 
In the present context an additional acceleration can be achieved by adjusting the 
accuracy of the approximation to the Neuberger-Dirac operator in the course of the 
inversion algorithm. This has first been described in ref. [51], and we now wish to 
present a second scheme that results in an even larger acceleration factor. 

We again consider eq. (8.1) in the chirality sector that does not contain the zero- 
modes. The adapted-precision algorithm that will be discussed should then be ap- 
plied to the preconditioned system (8.7), but in order to keep the presentation as 
simple as possible we shall stick to the original system in this section. 

9.1 Conjugate gradient algorithm 

For the basic inversion algorithm we choose the standard conjugate gradient method. 
To fix our notations we now first write down the defining equations of this algorithm 
even though these are very well known [28,45]. 

The algorithm generates a sequence tpi-, i = 1)2,3,..., of approximate solutions 
recursively together with an accompanying sequence pi of so-called search directions. 
If we introduce the residues 

g i = r j -Ai) i , (9.1) 
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the recursion is given by 

V>i+i =ipi + aip u a,i = {gi,9i)/{Pi,Api), (9.2) 
9i+i =9i~ aiApi, (9.3) 

Pi+i = 9i+i + biPi, h = (gi+i,gi+i)/(gi,9i), (9.4) 
and is usually started by setting 

^i = 0, gi = V, Pi = 9i- (9-5) 

Further information (on orthogonality and convergence properties in particular) can 
be found in the books quoted above. 

9.2 Iterative improvement 

Once an approximate solution xi °f the linear system is found, by performing a 
certain number of conjugate gradient iterations for example, it can be improved in 
the following way. First note that the norm of the residue 

771 = 77 - Axi (9.6) 

is much smaller than 1 1 77 1 1 as otherwise we would not consider \\ to be an approximate 
solution. Secondly it is obvious that the solution \ of the subtracted system 

A X = vi (9-7) 

yields the solution of the original system through 

i> = X i+X- (9-8) 

The important point to note is that x wm i n general be a small correction to xi since 
771 is small. In particular, if we solve the subtracted system to a relative precision 
of a few decimal places, the accuracy of the total field tp is improved by these many 
digits. 

We can now iterate this procedure and obtain the solution in the form 

00 

V> = J>fc, (9-9) 

fc=i 

Vo = V, Vk = Vk-i ~ Axk (A; = 1,2,3,...), (9.10) 
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where the fields Xk are assumed to be approximate solutions of the subtracted linear 
systems A\ = r/k-i- Depending on how accurate the solutions are, their magnitude 
llXfell will decrease more or less rapidly, but the algorithm is in any case always exact 
since the inaccuracies are eventually corrected for. 

9.3 Computation of \k 

So far no acceleration has been achieved and we have merely reorganized the calcu- 
lation in a fairly complicated way. The idea is now to use a numerical representation 
of A with low precision when calculating Xk- in doing so we should of course be very 
careful to choose a sufficiently accurate approximation as otherwise it is possible that 
the conjugate gradient algorithm becomes inefficient or even unstable. 

Before addressing this question, we recall that A acts in the subspace of negative 
chirality where 

2 

A = - P-D m ,P-, rri = \am 2 . (9.11) 

So if we use this representation of the operator and if a minmax polynomial is chosen 
that achieves a certain relative precision 5, the total approximation error is 

\\A-A\\ <t = + -\{am) 2 }{2{K + + K-)+5}. (9.12) 

Note incidentally that the original form (8.1) of the operator offers no numerical 
advantage at this point. The representation (9.11) is in fact more efficient by nearly 
a factor of 2 for a given uniform approximation error which is what counts here. 

Now when the conjugate gradient algorithm is applied to solve the linear system 
Ax = Vk-i, a sequence Xk,i of approximate solutions is generated (cf. subsect. 9.1). 
The decrease of the associated residues 



9k,i = Vk-i ~ Axk,i (9.13) 

can then be observed and the recursion is stopped as soon as ||<7fc,i|| falls below a 
certain level or if the number of iterations reaches a specified maximal value. In the 
course of this calculation a minimal requirement is that the residues are obtained to 
some precision at least, i.e. the chosen approximation for the operator A should be 
such that the error r||xfc,i|l is always smaller than ||<?fc,i|l- 

Since we cannot know in advance how large the approximate solutions \k,i are 
going to be, the accuracy r should be adjusted dynamically by increasing the degree 
of the minmax polynomial when needed. A relatively low precision, say 5 = 10 -3 , is 
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certainly good enough for the first few iterations, but as soon as r||x/c,i|| approaches 
||</fc,i||, the program should switch to a better approximation. In general this may 
happen several times before the recursion is halted and the calculated approximate 
solution Xk is returned to the calling program. 

9.4 Computation of rjk and global stopping criterion 

The residues rjk can in principle be calculated via eq. (9.10), but it is safer to obtain 
them directly from the current total solution vector through 



In this way r]k is guaranteed to be the true residue of the current approximation to 
the solution t/j of the linear system that we wish to solve eventually. In particular, 
if the global stopping criterion is based on the norm of the so computed residue, we 
need not worry about a possible accumulation of approximation errors in the course 
of the calculation. The algorithm simply stops when and only when an algorithm- 
independent criterion is fulfilled. 

It is not difficult to figure out how good the numerical approximation of A must be 
in eq. (9.14). Basically the approximation error r times the norm of the current total 
solution vector should be small compared to \\r]k\\ so that the residue is correctly 
obtained to a few decimal places. The stopping criterion for the next subtracted 
system, Ax = rjk, then has to be such that the accuracy of the solution is not driven 
to a point, where the random digits in the source rjk become important. 

9. 5 Fine tuning of the algorithm 

To achieve the maximal possible acceleration, the degree of the minmax polynomial 
should never become very large in the course of the solution of the subtracted sys- 
tems. It is therefore important that the conjugate gradient algorithm stops before 
this happens, which will be the case if the limits on the residue and the maximal 
number of iterations are properly chosen. 

A small additional improvement is obtained if the search direction is passed from 
the current subtracted system to the next as if the whole algorithm was a single 
sequence of conjugate gradient iterations. This is a correct procedure up to rounding 
and approximation errors that can partly be corrected by adjusting the initial search 
directions pk,i such that the orthogonality relations 



A- 




(9.14) 



1=1 



{(Pk,i ~ 9k,i),9k,i) = ((Pk,i ~ 9k,i),Ap k: i) = 



(9.15) 
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are restored to machine precision. 

If the stopping criteria are well chosen, this algorithm in general leads to significant 
savings in computer time (factors of 3 to 4 are not uncommon) . As in the case of the 
calculation of the index of the Dirac operator discussed in sect. 6, the use of 32 bit 
arithmetic can give another factor of 2 or so on current PC processors [48]. Single- 
precision arithmetic is in fact completely adequate for the approximate solution of 
the subtracted systems [but will usually be insufficient for the accumulation of the 
total solution vector and the calculation of the residues (9.14)]. 



10. Concluding remarks 

Numerical simulations in the e-regime of QCD are difficult, not only because they 
tend to require large amounts of computer time but also because special care needs 
to be taken to remain in control of the approximation errors. In this paper we have 
presented an algorithmic framework that is fully satisfactory in this respect and that 
makes the calculations significantly faster. 

Although we focused on a particular formulation of lattice QCD, it is clear that 
many of the methods that we have discussed are more widely applicable. An obvious 
case in this connection are the variants of lattice Dirac operators and gauge field 
actions that have recently been studied in refs. [18,52-56]. It is also conceivable 
that, under certain conditions, even ordinary formulations of lattice QCD can profit 
from mixed-precision algorithms and low-mode preconditioning. 

A disadvantage of our algorithms is perhaps that the propagators for different 
quark masses cannot be obtained simultaneously. In the e-regime the quark masses 
are very small, however, and apart from the zero-mode contribution (which must be 
calculated only once), the quark propagator should not vary a lot from one value of 
the quark mass to another. By taking linear combinations of previously computed 
propagators as a first approximation, the computer time needed for the calculation 
of the propagator at the next quark mass should thus be significantly reduced. 

The computation of the quark propagator in the e-regime, using just the standard 
conjugate gradient algorithm, requires a computational effort that grows roughly like 
the square of the lattice volume at large volumes. Once low-mode preconditioning 
is switched on, and if the number of subtracted modes is scaled proportionally to 
some fractional power of the volume, the situation becomes less transparent, but it 
is possible that a more favourable asymptotic behaviour will result. In any case the 
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total acceleration factor achieved by the methods described in this paper is quite 
impressive, and we hope that physically interesting computations in the e-regime 
will now be feasible on much larger lattices than previously considered. 

We are indebted to Pilar Hernandez, Mikko Laine, Laurent Lellouch and Peter 
Weisz for helpful discussions on chiral perturbation theory and on the distribution 
of the low-lying eigenvalues of the Dirac operator. We also wish to thank DESY for 
computer time and the staff of the computer centre for their support. L.G. was sup- 
ported in part by the EU under contract HPRN-CT-2000-00145. C.H. acknowledges 
support from the EU under contracts HPMF-CT-2001-01468 and HPRN-CT-2000- 
00145. 



Appendix A. Notations 

A.l Index conventions and Dirac matrices 

Lorentz indices are taken from the middle of the Greek alphabet and run from to 
3. Where appropriate repeated indices are automatically summed over. The Dirac 
matrices are assumed to satisfy 



but are otherwise left unspecified. It is often advantageous, however, to choose a 
chiral representation where 



In particular, pairs of Weyl fields (such as and in sect. 8) can then conveniently 
be stored in the memory space allocated for a Dirac field. 

A. 2 Wilson-Dirac operator 

In terms of the gauge-covariant forward and backward difference operators 



(7.) 



{7m , lu] = 2<5, 



(A.l) 




(A.2) 



V^{x) 




(A.3) 



V*V(*) 



1 



{tp(x) — U(x — a/t, fi) 1 ip(x — a/t)} 



(A.4) 



a 



26 



the Wilson-Dirac operator is given by 

^w = H7^ + V fl )-aVX}, (A.5) 

where a denotes the lattice spacing, U(x, fi) G SU(3) the link variables and ft the 
unit vector in direction fi. 

A. 3 Chebyshev polynomials 

The Chebyshev polynomials T n (z), n = 0, 1, 2, . . ., are defined through 

T„(cos0) = cos(n6>). (A.6) 

In particular, Tq{z) = 1 and T\(z) = z while for larger degrees the polynomials may 
be obtained algebraically using 

T n+l {z)-2zT n {z)+T n _ l {z) = 0, n>l. (A.7) 

The Clenshaw recursion [45] is a direct consequence of these equations. 



Appendix B. Proof of lemma 5.1 

We first note that 

||(1-P)QP|| = ||PQ(1-P)|| = 0. (B.l) 

This is a direct consequence of the definition (5.2) of g and of the hermiticity of the 
operators Q and P. If we define the projected operator 

Q = PQP+(1-P)Q(1-P), (B.2) 

it follows that 

\\Q - Q\\ = ||(1 - F)QF + PQ(1 -F)\\ = g (B.3) 

since P and 1 — P project to orthogonal subspaces. The min-max principle (ref. [49], 
Theorem XIII. 1) now implies that the eigenvalues of Q and Q, ordered in ascending 
order and counting multiplicities, are by no more than a distance g apart. 
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Appendix C. Proof of lemma 5.2 



We focus on the case of P + since the argument is exactly the same for the other 
projector. So let us consider an eigenvalue v k > and let u k £ V be the associated 
normalized eigenvector of PQP [eq. (5.3)]. From the definition (5.5) of g k we then 
infer that 

{Q-v k )u k = w k , \\w k \\ = Q k . (C.l) 
Next we write down the equations 

|| [1 - (P+) eX act] U k \\ 2 = (u k , [1 - (P+ ) exact ] U k ) 

= ((Q - AW, 1 ^ F _ + A)2 act (Q ~ A H) (C-2) 

that are trivially true as long as A is not in the spectrum of Q. We can actually 
safely set A = v k because the distance d k of v k from the spectrum in the range of the 
projector 1 — (P + ) exact does not vanish. The norm of the operator ratio in eq. (C.2) 
is then at most l/d 2 k and this leads to the inequality 

||[l-(P+)exact]u fe || < Q k /d k (C.3) 

when eq. (C.l) is taken into account. 

Together with the representation (5.4) of the projector P + in terms of the eigen- 
vectors u k , this bound implies 

||[l-(P+)exact]P+l>||< J2(e k /d k )\(u k ,v)\ (C.4) 

for any fermion field v, from which we deduce that 

||P+ [1 - (P+)exact] || = || [1 " (P+)exact] P+|| < «+■ (C.5) 

This is not quite what the lemma says, but the inequality is in fact more restrictive 
than it seems, because (P + ) exact and P + are, by definition, projectors of the same 
rank. In particular, the trace of the difference EI = P_|_ — (P+) exa ct vanishes, and a 
few lines of algebra show that 

M 2 +1HI = 2P + [1- (P+) e xact]P+ 

+ P+ [1 - (P+)exact] (1 " P+) + (1 " P+) [1 " (P+)exact] P+ • (C.6) 
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Fig. 2. Plot of f(h) — h 2 + h and graphical determination of the ranges of h (hashed 
bars) where the bounds (C.8) are satisfied. 



The operator thus satisfies 
H f =H, Tr H = 0, 



+ M|| < k + (1 + 2k+) 



(C.7) 



and consequently has to be of order k+. 

To make this explicit, we set k = k + (1 + 2k + ) and note that the eigenvalues h of 
H have to be such that 

-K<h 2 + h<K (C.8) 
(see fig. 2). In terms of k the condition stated on the last line of the lemma reads 

2(Z + 1)k<1 (C.9) 
so that k< \ and 



\h\< 



1 - 2k 



or \h+ 1| < 



K 



1-2k 



(CIO) 



Now since H acts in a subspace of dimension 21 at most, and since it has vanishing 
trace, the second range in (C.10) is excluded if 



(21 - 1) 



1 - 2k 



< 1 - 



1 - 2k 



(C.11) 
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which is the case in view of (C.9). All eigenvalues of EI thus satisfy the first inequality 
in (CIO) so that 

||H|| < T -^-. (C.12) 
Recalling the definition of M and k, this proves the lemma. 
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